UlC 

The  University  of  Illinois  at  Chicago 


A  QR  DECOMPOSITION  METHOD 
FOR  FLEXIBLE  MULTIBODY 
DYNAMICS 


Department  of  Mechanical  Engineering 


Approved  for  public  release; 
Diatributiop  UnBniil^ 


jfnc  QXJALrrr  inspected  a 


Technical  Report  #  MBS97-2-UIC 
Department  of  Mechanical  Engineering 
University  of  Illinois  at  Chicago 

November  1997 


A  QR  DECOMPOSITION  METHOD 
FOR  FLEXIBLE  MULTIBODY 
DYNAMICS 


19971215  083 


A.A.  Shabana 

Department  of  Mechanical  Engineering 
University  of  Illinois  at  Chicago 
842  West  Taylor  Street 
Chicago,  Illinois  60607-7022 


This  research  was  supported  by  the  U.S.  Army  Research  Office,  Research  Triangle  Park,  NC 


quality  IifOTf-frrlf!)  | 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  NO.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comment  regarding  this  burden  estimates  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services.  Directorate  for  information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington.  DC  20503. 


1 .  AGENCY  USE  ONLY  (Leave  blank) 


A.  TITLE  AND  SUBTITLE 


2.  REPORT  DATE 

November  1997 


3.  REPORT  TYPE  AND  DATES  COVERED 

Technical 


5.  FUNDING  NUMBERS 


A  QR  Decomposition  Method  for  Flexible  Multibody  Dynamics 


6.  AUTHOR(S) 


A.  A.  Shabana 


DAAG55-97-1-0303 


7.  PERFORMING  ORGANIZATION  NAMES(S)  AND  ADDRESS(ES) 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


University  of  Illinois  at  Chicago 
Chicago,  Illinois  60607-7022 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

U.S.  Army  Research  Office 
P.O.  Box  12211 

Research  Triangle  Park,  NC  27709-2211 


10.  SPONSORING  /  MONITORING 
AGENCY  REPORT  NUMBER 


AR0  35711.5-EG 


11.  SUPPLEMENTARY  NOTES 


The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and  should  not  be  construed  as 
an  official  Department  of  the  Army  position,  policy  or  decision,  unless  so  designated  by  other  documentation. 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 


12  b.  DISTRIBUTION  CODE 


Approved  for  public  release;  distribution  unlimited. 


13.  ABSTRACT  (Maximum  200  words) 


ABSTRACT  IN  TECHNICAL  REPORT 


15.  NUMBER  IF  PAGES 


16.  PRICE  CODE 


17  SECURITY  CLASSIFICATION  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRACT 
OR  REPORT  OF  THIS  PAGE  OF  ABSTRACT 


UNCLASSIFIED 


NSN  7540-01-280-5500 


UNCLASSIFIED 


UNCLASSIFIED 


standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  239-18 
298-102 


ABSTRACT 


Deformable  components  in  multibody  systems  are  subject  to  kinematic  constraints  the  represent 
mechanical  joints  and  specified  motion  trajectories.  These  constraints  can,  in  general,  be  described 
using  a  set  of  nonlinear  algebraic  equations  that  depend  on  the  system  generalized  coordinates  and 
time.  This  paper  describes  an  efficient  procedure  for  the  computer  implementation  of  the  absolute 
nodal  coordinate  formulation  for  flexible  multibody  applications.  In  the  absolute  nodal  coordinate 
formulation,  no  infinitesimal  or  finite  rotations  are  used  as  nodal  coordinates.  The  configuration  of 
the  finite  element  is  defined  using  global  displacement  coordinates  and  slopes.  By  using  this  mixed 
set  of  coordinates,  beam  and  plate  elements  can  be  treated  as  isoparametric  elements.  As  a 
consequence,  the  dynamic  formulation  of  these  widely  used  elements  using  the  absolute  nodal 
coordinate  formulation  leads  to  a  constant  mass  matrix.  It  is  the  objective  of  this  study  to  develop  a 
an  efficient  computational  procedure  that  exploits  this  feature.  In  this  procedure,  an  optimum  sparse 
matrix  structure  is  obtained  for  the  deformable  bodies  using  the  QR  decomposition.  Using  the  fact 
that  the  element  mass  matrix  is  constant,  a  QR  decomposition  of  a  modified  constant  connectivity 
Jacobian  matrix  is  obtained  for  the  deformable  body.  A  constant  velocity  transformation  is  used  to 
obtain  an  identity  generalized  inertia  matrix  associated  with  the  second  derivatives  of  the  coordinates 
used  in  the  absolute  nodal  coordinate  formulation,  thereby  minimizing  the  number  of  non-zero  entries 
of  the  coefficient  matrix  that  appears  in  the  augmented  Lagrangian  formulation  of  the  equations  of 
motion  of  the  flexible  multibody  systems.  The  computational  procedure  proposed  in  this  investigation 
can  be  used  for  the  treatment  of  large  deformation  problems  in  flexible  multibody  systems.  It  has  also 
the  advantages  of  the  algorithms  based  on  the  floating  frame  of  reference  formulations  since  it  allows 
for  easy  addition  of  general  nonlinear  constraint  and  force  functions. 


KEY  WORDS: 


Multibody  Dynamics,  Finite  Element  Method,  QR  Decomposition,  Absolute 
Nodal  Coordinate  Formulation,  Floating  Frame  of  Reference  Formulation, 
Incremental  Methods,  Large  Deformation,  Large  Rotation. 


1. 


INTRODUCTION 


The  performance  and  efficiency  of  multibody  simulation  codes  depend  largely  on  the  selection  of  the 
coordinates  used  to  formulate  the  dynamic  equations  of  multibody  systems  (Kim  and  Vanderpoleg, 
1986;  Mani  et  al.,  1985;  and  Singh  and  Likins,  1985).  The  choice  of  the  coordinates  defines  the 
structure  of  the  system  equations  of  motion  as  well  as  the  numerical  procedure  required  for  the 
solution  of  these  equations.  Extensive  research  efibrts  have  been  devoted  to  examine  the  effect  of  the 
coordinate  selection  on  the  complexity  of  the  formulation  as  well  as  the  efficiency  and  performance 
of  the  computer  algorithms  used  in  rigid  multibody  dynamics.  The  research  in  flexible  multibody 
dynamics,  on  the  other  hand,  has  been  primarily  focused  on  some  fundamental  issues  related  to 
modeling  the  dynamic  motion  of  flexible  bodies  that  undergo  large  displacements  (Shabana,  1997; 
and  Wasfy  and  Noor,  1997).  Several  finite  element  formulations  have  been  proposed  for  the  large 
displacement  analysis  of  flexible  multibody  systems.  Among  these  formulations  are  the  floating  firame 
of  reference  method  (Shabana,  1989),  the  incremental  methods  (Bel3^schko  and  Hsieh,  1973;  Rankin 
and  Brogan,  1986),  and  large  rotation  vector  formulations  (Simo  and  Vu-Quoc,  1986).  In  the  floating 
frame  of  reference  formulation,  which  is  the  most  widely  used  method  for  flexible  multibody 
dynamics,  two  sets  of  coordinates  are  used  to  define  the  configuration  of  the  flexible  body.  The  first 
set  is  the  set  of  reference  coordinates  which  defines  the  location  and  orientation  of  a  selected 
deformable  body  coordinate  system.  The  second  set  is  the  set  of  elastic  coordinates  which  defines  the 
deformation  of  the  body  with  respect  to  its  coordinate  system.  The  floating  firame  of  reference 
formulation  leads  to  a  highly  nonlinear  mass  matrix  as  the  result  of  the  inertia  coupling  between  the 
rigid  body  motion  and  the  elastic  deformation.  This  formulation,  however,  can  be  used  to  obtain  an 
exact  representation  of  the  rigid  body  dynamics  and  leads  to  zero  strains  under  an  arbitrary  rigid  body 
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motion,  even  in  the  case  when  non-isoparametric  finite  elements  such  as  conventional  beam  and  plate 
elements  are  used. 

Incremental  finite  element  formulations  have  been  successfully  used  in  the  large  deformation 
analysis  of  structural  systems.  In  the  incremental  methods,  the  configuration  of  the  finite  element  is 
described  using  the  element  nodal  coordinates.  Since  non-isoparametric  elements  can  not  describe 
an  arbitrary  large  rotation,  the  large  rotation  of  the  element  is  represented  as  sequence  of  infinitesimal 
rotations.  The  infinitesimal  rotations  can  then  be  described  using  the  conventional  element  shape 
function  and  the  element  nodal  coordinates.  It  is  important  to  note,  however,  that  when  non- 
isoparametric  elements  such  as  beams  and  plates  are  used,  the  incremental  methods  can  not  be  used 
to  obtain  exact  representation  of  the  rigid  body  dynamics,  and  such  methods  do  not  lead  to  zero 
strains  under  an  arbitrary  rigid  body  motion.  For  this  reason,  the  incremental  methods  are  not  widely 
used  in  flexible  multibody  computer  programs. 

The  large  rotation  vector  formulations  are  non-incremental  and  were  developed  to  circumvent 
the  partial  hnearization  used  in  the  finite  element  incremental  formulations.  In  the  large  rotation 
vector  formulations,  absolute  coordinates  and  finite  rotations  are  used  to  define  the  element 
configuration.  These  formulations  lead  to  a  simpler  inertia  matrix  and  a  more  complex  stiflfiiess 
matrix  Nonetheless,  the  mass  matrix  remains  nonlinear  when  three  dimensional  elements  are  used. 
Large  rotation  vector  formulations  also  lead  to  excessive  shear  forces  as  the  result  of  the  description 
used  for  the  finite  rotation  of  the  element  cross  section  (Shabana,  1997). 

As  previously  pointed,  the  floating  fi-ame  of  reference  formulation  is  the  most  widely  used 
method  in  the  dynamic  analysis  of  flexible  multibody  systems.  This  formulation  has  been  implemented 
in  several  commercial  and  research  computer  programs.  The  popularity  of  this  method  is  atributed 
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to  the  fact  that  general  constraints  and  forcing  functions  can  be  systematically  and  easily  introduced 
to  the  formulation  and  can  also  be  implemented  in  the  computer  programs  in  a  straightforward 
manner.  However,  the  use  of  the  floating  frame  of  reference  formulation  has  been  limited  to  small 
deformation  problems;  a  limitation  that  arises  as  the  result  of  the  t5?pe  of  coordinates  and  the  motion 
description  used  in  this  formulation.  Recently,  a  new  procedure  called  the  absolute  nodal  coordinate 
formulation  was  introduced  (Shabana,  1997).  In  this  formulation,  a  new  set  of  finite  element 
coordinates  is  employed.  This  set  of  coordinates  consists  of  global  displacement  coordinates  and 
slopes.  Using  this  set  of  coordinates,  beam  and  plate  elements  can  be  treated  as  isoparametric 
elements,  and  as  a  consequence,  exact  modeling  of  the  rigid  body  dynamics  can  be  obtained  using  the 
element  shape  function  and  the  element  nodal  coordinates.  Furthermore,  the  absolute  nodal 
coordinate  formulation  leads  to  zero  strains  under  an  arbitrary  rigid  body  motion.  When  such  a 
formulation  is  used,  some  of  the  fundamental  problems  encountered  when  using  the  floating  frame 
of  reference  formulation  can  be  avoided.  One  of  these  problems  is  the  selection  of  the  deformable 
body  coordinate  system  (Shabana,  1989).  This  problem  does  not  arise  in  the  absolute  nodal 
coordinate  formulation  since  the  coordinates  used  in  this  formulation  are  all  defined  in  the  global 
system.  The  absolute  nodal  coordinate  formulation  also  leads  to  a  constant  mass  matrix,  and  as  such, 
an  efficient  procedure  for  implementing  this  formulation  in  flexible  multibody  computer  programs  can 
be  developed. 

In  the  general  computational  algorithms  developed  for  the  nonlinear  dynamic  analysis  of 
multibody  systems,  nonlinear  constraint  equations  that  describe  mechanical  joints  and  specified 
motion  trajectories  can  be  systematically  introduced  to  the  dynamic  formulation.  This  provides  the 
flexibility  of  building  computer  models  for  mechanical  systems  with  complex  topological  structure. 
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By  utilizing  sparse  matrix  algebra,  these  models  can  be  efficiently  simulated  and  modified.  In  order 
to  maintain  the  generality  of  the  djmamic  formulation,  several  of  the  general  purpose  computer 
algorithms  developed  for  the  simulation  of  multibody  applications  are  based  on  solving  the  following 
system  of  equations  (Shabana,  1994); 


M 

q 

q; 

0 

X 

(1) 


where  M  is  the  mass  matrix  of  the  system,  is  the  Jacobian  matrix  of  the  kinematic  constraints,  q 
is  the  vector  of  the  system  generalized  coordinates,  X  is  the  vector  of  Lagrange  multipliers,  is  the 
vector  offerees  that  include  external,  gravity,  Coriolis,  centrifiigal,  and  elastic  forces,  and  is  the 
vector  resulting  from  the  differentiation  of  the  constraint  equations  twice  with  respect  to  time.  For 
instance,  the  constraint  equations  can  be  written  in  the  following  form: 

C(q,0  =  0  (2) 

where  C  is  the  vector  of  constraint  functions,  and  t  is  time.  Upon  differentiating  the  constarint 
equations  twice  with  respect  to  time,  one  obtains 

C,q=Q.  (3) 


where  is  the  vector  defined  as 

Q.=  -C„-(C,q)^q-2C„4  (4) 

The  coefficient  matrix  of  Eq.  1  has  a  sparse  matrix  structure.  An  efficient  solution  of  this  system  of 
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equations  can  be  obtained  if  the  number  of  non-zero  entries  of  the  nonlinear  coefficient  matrix  is 
minimized. 

It  is  the  objective  of  this  investigation  to  develop  a  new  computational  procedure  for  flexible 
multibody  dynamics  that  exploits  the  sparse  matrix  structure  presented  in  Eq.  1 .  The  new  procedure 
is  based  on  the  absolute  nodal  coordinate  formulation  which  leads  to  a  constant  mass  matrix  for  the 
deformable  body  (Shabana,  1997).  In  the  absolute  nodal  coordinate  formulation,  global  displacement 
coordinates  and  slopes  are  used  to  define  the  element  configuration.  Connectivity  conditions  between 
the  finite  elements  used  in  the  discretization  of  the  deformable  body  can  be  described  using  a  set  of 
linear  algebraic  constraint  equations.  Using  the  constant  Jacobian  matrix  of  these  constraint  equations 
and  the  fact  that  the  system  mass  matrix  is  constant,  an  efficient  QR  decomposition  for  a  constant 
generalized  Jacobian  matrix  can  be  obtained  only  once  before  the  dynamic  simulation  (Atkinson, 
1987;  Press  et  al.,  1992;  and  Strang,  1988).  Orthogonal  vectors  resulting  fi-omthis  decomposition 
are  used  to  define  an  identity  inertia  matrix  associated  with  the  coordinates  of  the  deformable  body, 
thereby  minimizing  the  number  of  non-zero  entries  of  the  sparse  coefficient  matrix  that  appears  in  Eq. 
1. 

In  the  analysis  presented  in  this  investigation  we  distinguish  between  two  sets  of  constraints. 
The  first  set  consists  of  the  constraints  that  describe  the  connectivity  between  the  elements  of  a 
deformable  body  in  the  multibody  system.  These  constraints  are  assumed  to  be  linear  functions  of  the 
nodal  coordinates  of  the  deformable  body.  We  will  refer  to  these  constraints  as  the  connectivity 
constraints.  The  forces  resulting  fiom  these  connectivity  constraints  will  be  systematically  eliminated 
from  the  equations  of  motion  of  the  deformable  body  using  the  QR  decomposition  of  the  a 
generalized  constraint  Jacobian  matrix.  The  second  set  of  constraints  consists  of  the  constraints  that 
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describe  the  joints  between  deformable  bodies  as  well  as  the  constraints  that  describe  specified  motion 
trajectories.  The  forces  of  these  constraints  will  not  be  eliminated  from  the  dynamic  formulation  in 
order  to  increase  the  generality  of  the  algorithm  presented  in  the  paper. 

This  paper  is  organized  as  follows.  In  Section  2,  the  absolute  nodal  coordinate  formulation 
is  briefly  reviewed  and  the  form  of  the  constant  inertia  matrix  associated  with  the  coordinates  used 
in  this  formulation  is  defined.  In  Section  3,  the  equations  of  motion  of  the  finite  elements  expressed 
in  terms  of  the  element  connectivity  forces  are  presented.  These  element  equations  are  used  to  define 
the  equations  of  motion  of  the  deformable  body.  In  Section  4,  the  constraint  equations  that  describe 
the  coimectivity  of  the  elements  of  the  deformable  body  are  presented  and  used  to  define  the 
connectivity  constraint  Jacobian  matrix.  Using  this  constraint  Jacobian  matrix,  the  element 
coimectivity  forces  are  expressed  in  terms  of  Lagrange  multipliers.  Utilizing  the  fact  that  the  mass 
matrix  of  the  deformable  body  is  constant,  a  generalized  coimectivity  Jacobian  matrix  is  defined  and 
its  QR  factors  are  presented  in  Section  5.  A  computational  procedure  for  eliminating  the  connectivity 
forces  from  the  equations  of  motion  of  the  deformable  body  is  presented  in  Section  6.,  while  the  final 
form  of  the  equations  of  motion  and  the  solution  algorithm  of  these  equations  are  presented  in 
Section  7.  Summary  and  discussion  of  the  work  described  in  this  paper  is  presented  m  Section  8. 

2.  ABSOLUTE  NODAL  COORDINATE  FORMULATION 

In  the  analysis  presented  in  this  paper,  two  dimensional  beam  elements  are  used  as  examples.  The 
procedure  developed  in  this  investigation,  however,  can  be  also  applied  to  three  dimensional  beam, 
plate  and  shell  elements.  Figure  1  shows  a  two-dimensional  beam  element  which  has  two  nodes 
defined  by  the  points  A  and  B.  In  the  absolute  nodal  coordinate  formulation,  global  displacements 
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coordinates  and  slopes  are  used  as  the  nodal  coordinates.  In  this  formulation,  no  infinitesimal  or  finite 
rotations  are  used  as  nodal  coordinates.  By  using  the  global  displacement  coordinates  and  slopes,  and 
a  shape  function  that  has  a  complete  set  of  rigid  body  modes;  beam  and  plate  elements  can  be  treated 
as  isoparametric  elements  that  can  be  used  to  obtain  exact  modeling  of  the  rigid  body  dynamics 
(Shabana  and  Christensen,  1997).  In  the  absolute  nodal  coordinate  formulation,  the  global  position 
of  an  arbitrary  point  on  an  element  j  on  the  deformable  body  i  in  the  multibody  system  can  be  written 
as 

i^  =  S^ei''  (5) 

where  is  the  global  position  vector  of  an  arbitrary  point  on  the  beam  element.  S'-'  is  the  element 
shape  function  which  has  a  complete  set  of  rigid  body  modes,  and  is  the  vector  of  absolute  nodal 
coordinates  and  slopes  of  the  element.  Using  the  preceding  equation,  the  kinetic  energy  of  the  finite 
element  can  be  defined  as 

yV 


where  and  are  the  mass  density  and  volume  of  the  finite  element  j  of  the  deformable  body  /. 
Using  the  preceding  two  equations,  it  can  be  shown  that  the  kinetic  energy  of  the  finite  element  can 
be  written  as 

(7) 

where  M'^  is  the  constant  mass  matrix  of  the  element  defined  as 
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(8) 


v« 


The  fact  that  the  element  mass  matrix  is  constant  plays  a  fundamental  role  in  the  algorithm  developed 
in  this  investigation.  Note  that  using  the  absolute  nodal  coordinate  formulation,  the  mass  matrix 
remains  also  constant  when  three  dimensional  beams  are  considered.  As  a  consequence,  the  vector 
of  Coriolis  and  centrifijgal  forces  is  identically  equal  to  zero.  The  stiffness  matrix,  on  the  other  hand, 
is  a  highly  nonlinear  function  of  the  element  coordinates,  even  in  the  case  when  small  deformation 
problems  are  considered. 

Since  in  the  absolute  nodal  coordinate  formulation,  the  coordinates  are  defined  in  the  global 
system;  there  is  no  reason  to  justify  using  different  interpolating  polynomials  for  the  displacement 
components.  For  example,  in  the  case  of  the  beam  element  shown  in  Fig.  1,  cubic  polynomials  can 
be  used  to  describe  both  components  of  the  displacements  defined  in  Eq.  5.  In  this  case,  the  element 
shape  function  can  be  defined  as 


S^'= 


1-3  ^+2^ 
0 


0  0 
i-se+ie  0  i{^-2e*e) 


3^-2^  0  Kf-e)  0 

0  3^-24®  0 


(9) 


where  /  is  the  length  of  the  element,  ^  =  x/l,  and  x  is  the  axial  coordinate  that  defines  the  position  of 
an  arbitrary  point  on  the  element  in  the  undeformed  state.  Using  this  shape  function,  the  vector  of 
nodal  coordinates  of  the  element  can  be  defined  as 
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kV-L'^  o'^  o'^  o'^ 

,  -jCj  ^2  *^3  M  *^5  ^6  ^1  *^i\ 


where  and  are  the  absolute  displacement  coordinates  of  the  node  at  A,  and  e/  are  the 
absolute  displacement  coordinates  of  the  node  at  B,  and 


3rj  (x=0) 


3/-2(x=0) 


8rj  (x=/) 

dx 


drJx=l) 


and  and  rj  are  the  components  of  the  vector  r.  Using  the  shape  function  of  Eq.  9,  the  constant  mass 
matrix  of  the  element  can  be  evaluated,  using  Eq.  8,  as 


M^'  -  m 


il 

0 

11/ 

0 

9 

0 

-  13/ 

0 

35 

210 

70 

“420 

0 

13 

0 

11/ 

0 

9 

0 

_  13/ 

35 

210 

70 

"*420 

11/ 

0 

/2 

0 

13/ 

0 

.  /' 

0 

210 

105 

420 

140 

0 

11/ 

0 

f 

0 

13/ 

0 

.  /' 

210 

105 

420 

140 

9 

0 

13/ 

0 

13 

0 

_  11/ 

0 

70 

420 

35 

210 

0 

9 

0 

13/ 

0 

13 

0 

_  11/ 

70 

420 

35 

210 

^  13/ 

0 

. 

0 

.  11^ 

0 

/2 

0 

420 

’l40 

^*210 

105 

0 

^  13/ 

0 

_  /2 

0 

,  11/ 

0 

P 

~420 

140 

”210 

105 

It  can  be  demonstrated  that  the  absolute  nodal  coordinate  formulation  leads  to  exact  modeling  of  rigid 
body  dynamics  and  does  not  lead  to  the  linearization  of  the  equations  of  motion  as  in  the  case  of 
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incremental  formulations.  The  equivalence  of  the  absolute  nodal  coordinate  formulation  and  the 
floating  frame  of  reference  formulation  can  be  demonstrated  if  local  slopes  instead  of  infinitesimal 
rotations  are  used  as  the  local  elastic  coordinates  in  the  floating  frame  of  reference  formulation.  This 
equivalence  was  recently  demonstrated  by  Shabana  and  Schwertassek  (1997). 

3.  FINITE  ELEMENT  EQUATIONS  OF  MOTION 

In  the  dynamic  formulation  presented  in  this  paper,  we  consider  two  types  of  constraints.  The  first 
type  of  constraints  is  the  constraints  that  represent  the  connectivity  conditions  between  the  finite 
elements  of  a  deformable  body.  These  constraints  are  linear  fiinctions  of  the  element  nodal 
coordinates.  The  second  type  of  constraints  is  the  nonlinear  constraints  that  represent  the  joints 
between  different  deformable  bodies,  or  specified  motion  trajectories.  These  nonlinear  constraints  are 
described  by  Eq.  2,  and  will  be  introduced  to  the  dynamic  formulation  using  the  technique  of 
Lagrange  multipliers  in  order  to  maintain  the  generality  of  the  formulation. 

The  first  type  of  linear  constraints  that  describe  the  connectivity  of  the  elements  of  one 
deformable  body  will  be  the  subject  of  discussion  in  this  and  the  following  sections.  A  non- 
conventional  algorithm  that  utilizes  the  fact  that  the  element  mass  matrix  is  constant  will  be  presented 

N 

and  its  advantages  will  be  discussed.  To  this  end,  the  equations  of  motion  of  the  finite  element  j  on 
the  deformable  body  i  are  written  as 

M^e'^'  =  F'i;+F'^  (13) 

where  F/  is  the  vector  of  element  forces  that  include  the  externally  applied  and  gravity  forces,  and 
the  elastic  forces  that  result  from  the  element  deformation,  and  F/  is  the  vector  of  constraint  forces 
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resulting  from  the  connectivity  of  the  element  with  the  other  elements  used  in  the  finite  element 
discretization  of  the  deformable  body.  Using  the  preceding  equation,  the  equations  of  motion  of  the 
deformable  body  i  can  be  written  in  terms  of  the  connectivity  forces  as 

M'e'  =  F;+F'  (14) 


where  M'  is  the  mass  matrix  of  the  deformable  body  given  in  terms  of  the  mass  matrices  of  its 
elements  as 


M’  = 


M'^  0 


M" 


(15) 


in  which  is  the  total  number  of  elements  of  the  deformable  body.  The  vectors  e',  F^',  and  F^'  are 
defined  as 


(16) 


Note  that  in  these  equations,  the  elements  are  not  assembled  and  none  of  the  element  coordinates  is 
eliminated.  The  equations  of  the  system,  however,  are  still  valid  since  the  connectivity  forces  appear 
in  these  equations. 
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4.  JACOBIAN  OF  THE  CONNECTIVITY  CONSTRAINTS 

The  constraints  that  describe  the  connectivity  between  the  finite  elements  of  the  deformable  body  i 
will  be  formulated  in  this  investigation  using  a  set  of  algebraic  equations  that  depend  linearly  on  the 
nodal  coordinates.  In  this  case,  these  constraint  equations  can  be  written  in  the  following  form; 

0'(e')  =  c  (17) 

where  O'  is  the  connectivity  constraint  vector,  and  c  is  a  constant  vector.  Note  that  the  connectivity 
constraints  do  not  depend  on  time.  For  a  virtual  change  of  the  body  coordinates,  the  preceding 
equation  leads  to 

O'6e'  =  0  (18) 

where  is  the  Jacobian  matrix  of  the  connectivity  constraints.  The  connectivity  constraint  forces 
in  Eq.  14  can  be  expressed  in  terms  of  the  Jacobian  matrix  of  the  element  connectivity  constraints  as 

F' x;  (19) 

where  XJ  is  the  vector  of  Lagrange  multipliers  associated  with  the  connectivity  constraints  of  the 
deformable  body.  Substituting  Eq.  19  into  Eq.  14,  one  obtains 

M'e'  =  F;-Of  (20) 

Since  the  mass  matrix  M'  is  constant  and  block-diagonal,  its  inverse  can  be  efficiently  calculated 
onlynce  before  the  dynamic  simulation  and  used  to  write  the  preceding  equation  in  the  following 
form; 
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(21) 


which  can  be  written  as 

e'  =  (M')*'Fl-A'A'  (22) 

where 

A'=(M')'^Of  (23) 

is  a  constant  generalized  Jacobian  matrix.  This  matrix,  which  needs  to  be  evaluated  only  once  before 
the  dynamic  simulation,  plays  a  significant  role  in  the  computational  procedure  proposed  in  this  paper. 


5.  FACTORS  OF  THE  GENERALIZED  JACOBIAN 

The  connectivity  constraint  forces  acting  on  the  elements  of  the  deformable  body  can  be 
systematically  eliminated  using  the  QR  factors  of  the  generalized  Jacobian  matrix  of  Eq.  23.  Using 
thq  QR  decomposition,  the  generalized  Jacobian  matrix  of  Eq.  23  can  be  written  as 

A'  =  QR  (24) 

where  Q  is  an  orthogonal  matrix  (Q^  Q  =  I ),  and  R  is  an  upper-triangular  matrix.  If  A'  is  an  w  x  n 
matrix,  then  Q  is  an  w  x  «  and  R  is  n  x  »  matrix.  Since  the  number  of  connectivity  constraints  is 
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always  smaller  than  the  number  of  coordinates,  m'^n.la  this  case,  the  preceding  equation  can  be 
written  as  (Kim  and  Vanderploeg,  1986) 


A'=[Q.  QJ 


R, 

0 


(25) 


where  Q,  and  Qj  are  the  partitions  of  the  orthogonal  matrix  Q.  The  dimension  of  Qj  is  /»  x  n,  while 
the  dimension  of  Q2  is  w  x  ( m  -  w ).  Since  Q  is  an  orthogonal  matrix,  it  follows  that 

QIQ,  =  0  (26) 

It  is  also  clear  from  Eq.  25  that 

A'  =  QjRi  (27) 


Using  the  preceding  three  equations,  it  can  be  verified  that 

Q2A'  =  0  (28) 

This  result  will  be  used  to  eliminate  the  connectivity  forces  from  the  equations  of  motion  of  the 
deformable  body  as  demonstrated  in  the  following  section. 


6.  ELIMINATION  OF  THE  CONNECTIVITY  FORCES 

Recall  that  the  matrix  A'  is  a  constant  matrix,  and  therefore,  the  matrix  Q2  that  results  from  the  QR 
decomposition  of  A!  is  also  constant.  This  fact  can  be  utilized  to  obtain  an  efficient  procedure  for 
eliminating  the  connectivity  forces  and  define  an  identity  generalized  inertia  matrix  for  the  deformable 
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body.  To  this  end,  the  following  velocity  transformation  is  employed: 

e'  =  Q2q'  (29) 

Substituting  this  coordinate  transformation  into  Eq.  18,  and  premultiplying  by  the  transpose  of  the 

matrix  Qj,  one  obtains 

QlQ2q'=Ql(M')'‘Fl-Q,A'Xi  (30) 

Since  the  columns  of  the  matrix  Q2  represent  a  set  of  orthogonal  vectors,  one  has 

qJq,  =  I  (31) 

Substituting  Eqs.  28  and  31  into  Eq.  30,  one  obtains 

q'  =  Qj(M')''Fl  (32) 

which  shows  that  the  generalized  inertia  matrix  associated  with  the  set  of  coordinates  q'  is  the  identity 
matrix.  As  a  consequence,  a  number  of  non-zero  entries  equal  to  the  number  of  the  coordinates  in  the 
set  q'  needs  to  be  stored  in  order  to  define  the  inertia  matrix  of  the  deformable  body.  Note  also  that 
in  Eq.  32,  the  forces  of  the  connections  between  the  finite  elements  of  the  deformable  body  are 
eliminated.  Equation  32,  therefore,  reduces  to  the  simple  form 

q'  =  Qi  (33) 

where  Qj  is  the  vector  of  generalized  forces  defined  as 
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Q;  =  Qj(M')‘‘Fi 


(34) 


Equation  33  is  in  a  form  suitable  for  implementation  in  general  purpose  flexible  multibody  computer 
programs  that  e^qjloit  sparse  matrix  techniques. 


7.  SYSTEM  EQUATIONS  AND  NUMERICAL  PROCEDURE 
Non-linear  constraints  between  deformable  bodies  can  be  formulated  in  the  absolute  nodal  coordinate 
formualation  using  Eq.  2.  In  this  case,  the  augmented  form  of  the  equations  of  motion  takes  the 
following  form: 


I  < 

q 

q; 

X 

(35) 


where  in  this  case  the  generalized  inertia  matrix  of  the  system  reduces  to  an  identity  matrix.  In  Eq. 
35,  C,  is  the  Jacobian  matrix  of  the  kinematic  constraints  that  describe  the  joints  between  deformable 
bodies  as  well  as  specified  motion  trajectories.  The  vector  q  in  Eq.  35  is  the  vector  of  the  system 
coordinates  which  can  be  written  as 


q  = 


,2T 


where  q'  is  the  vector  of  generalized  coordinates  of  the  deformable  body  i,  and  is  the  total  number 
of  bodies  in  the  system.  Equation  35  can  be  efficiently  solved  for  the  accelerations  and  Lagrange 
multipliers  using  sparse  matrix  techniques  (Duff  et  al.,  1986;  Press  et  al.,  1992).  The  vector  of 
Lagrange  multiphers  can  be  used  to  determine  the  generalized  joint  forces.  Using  the  constraint 
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Jacobian  matrix  of  the  joint  constraints  C, ,  a  set  of  independent  coordinates  can  be  identified.  The 
accelerations  associated  with  the  independent  coordinates  can  be  integrated  forward  in  time  in  order 
to  determine  the  independent  velocities  and  coordinates.  Using  the  independent  coordinates,  Eq.  2 
can  be  solved  for  the  dependent  coordinates  using  a  Newton-Raphson  algorithm.  Dependent 
velocities  can  be  determined  using  the  velocity  kinematic  relationships  (Shabana,  1994) 

8.  SUMMARY  AND  DISCUSSION 

Several  finite  element  formulations  have  been  proposed  for  the  dynamic  analysis  and  simulation  of 
flexible  multibody  systems.  Among  these  formulations  are  the  floating  frame  of  reference  method,  the 
incremental  methods,  and  the  large  rotation  vector  formulations.  The  floating  frame  of  reference 
formulation  is  the  most  widely  used  method  for  the  dynamic  analysis  of  flexible  multibody  systems 
since  such  a  formulation  allows  easy  addition  of  general  constraint  and  force  functions.  The  use  of 
this  method,  however,  has  been  limited  to  small  deformation  problems  because  of  the  nature  of  the 
generalized  coordinates  used.  Incremental  methods,  on  the  other  hand,  had  less  acceptance  in  the 
multibody  community  because  of  the  linearized  kinematic  equations  used  to  describe  the  overall 
motion  of  the  finite  element  when  these  methods  are  used.  Because  of  this  kinematic  description, 
exact  modding  of  the  rigid  body  dynamics  can  not  be  obtained  using  the  incremental  methods  when 
non-isoparametric  finite  elements  such  as  conventional  beam  and  plate  elements  are  used.  The 
coordinates  used  in  the  non-incremental  large  rotation  vector  formulations  include  finite  rotations 
which  are  described  using  interpolation  polynomials  in  a  similar  manner  to  the  displacement 
coordinates.  Such  a  motion  description  leads  to  excessive  shear  forces  that  lead  to  serious  numerical 
and  fundamental  modeling  problems. 
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In  this  investigation,  a  new  computational  finite  element  procedure  is  developed  for  the 
computer-aided  analysis  of  flexible  multibody  systems.  This  procedure,  which  is  based  on  the  absolute 
nodal  coordinate  formulation,  leads  to  an  optimum  sparse  matrix  structure  and  allows  for  easy 
addition  of  constraints  and  forcing  functions,  thereby  maintaining  the  main  advantages  of  the 
algorithms  based  on  the  floating  frame  of  reference  formulation.  Furthermore,  the  new  procedure  can 
be  used  for  the  large  deformation  analysis  of  flexible  multibody  systems,  and  as  such,  it  does  not 
suffer  fi-om  the  limitation  of  the  floating  firame  of  reference  formulation.  In  the  absolute  nodal 
coordinate  formulation,  global  displacement  coordinates  and  slopes  are  used  to  define  the  element 
configuration.  Infinitesimal  or  finite  rotations  are  not  used  as  nodal  coordinates.  By  using  this  set  of 
coordinates,  exact  modeling  of  the  rigid  body  dynamics  can  be  obtained  using  the  element  shape 
function  and  the  nodal  coordinates.  As  a  consequence  of  using  these  coordinates,  beam  elements  can 
be  treated  as  isoparametric  elements  and  an  arbitrary  rigid  body  displacement  leads  to  zero  strains. 
The  absolute  nodal  coordinate  formulation  also  leads  to  a  constant  mass  matrix,  and  as  a  result,  the 
vector  of  Coriolis  and  centrifugal  forces  is  identically  equal  to  zero. 

In  this  investigation,  an  efficient  implementation  of  the  absolute  nodal  coordinate  formulations 
for  flexible  multibody  applications  is  described.  In  the  procedure  described  in  this  paper,  advantage 
is  taken  of  the  fact  that  the  element  inertia  matrix  is  constant.  Two  different  sets  of  kinematic 
constraints  are  considered.  The  first  set  consists  of  the  constraints  that  describe  the  connectivity  of 
the  elements  of  a  deformable  body  in  the  multibody  system.  These  constraints  which  are  linear 
functions  of  the  nodal  coordinates  are  systematically  eliminated  using  the  procedure  described  in  this 
paper.  The  second  set  consists  of  the  constraints  that  describe  mechanical  joints  between  deformable 
bodies  in  the  system  as  well  as  specified  motion  trajectories.  These  constraints  are,  in  general. 
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nonlinear  and  they  are  not  eliminated  from  the  dynamic  formulation  in  order  to  increase  the  generality 
of  the  procedure  described  in  this  paper.  The  equations  of  motion  of  the  finite  elements  are  first 
formulated  in  terms  of  a  redundant  set  of  coordinates.  These  equations  are  expressed  in  terms  of  the 
connectivity  forces  which  are  expressed  in  terms  of  the  constraint  Jacobian  matrix  and  lagrange 
multipliers.  A  constant  generalized  Jacobian  matrix  is  obtained  using  the  inverse  of  the  constant 
element  mass  matrix.  A  QR  decomposition  is  used  to  determine  the  factors  of  the  generalized 
Jacobian  matrix.  Using  the  orthogonal  matrix  in  the  QR  decomposition,  a  velocity  transformation 
which  consists  of  orthogonal  column  vectors  is  developed  and  used  to  express  the  coordinates  of  the 
elements  of  the  deformable  body  in  terms  of  a  reduced  set  of  generalized  coordinates.  The  inertia 
matrix  associated  with  the  reduced  set  of  coordinates  is  the  identity  matrix.  Using  this  result,  an 
optimum  sparse  matrix  structure  for  the  equations  of  motion  of  the  flexible  multibody  system  can  be 
defined.  The  resulting  system  of  equations  of  motion  can  be  solved  using  the  numerical  procedure 
described  in  Section  7.  Note  also  that  despite  the  fact  that  the  element  connectivity  constraint  forces 
are  eliminated  from  the  final  form  of  the  equation  of  motion,  these  forces  can  be  easily  computed 
using  the  procedure  described  in  this  paper.  To  this  end,  the  nodal  accelerations  of  the  deformable 
body  i  can  be  computed  from  the  generalized  accelerations  of  the  body  using  the  time  derivative  of 
Eq.  29.  These  nodal  accelerations  can  be  substituted  in  Eq.  14  in  order  to  determine  the  connectivity 
force  vector  F/. 
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